#include "../includes/need-a-name.hh"
#include "includes/mhd.hh"
#include <map>

MHD::MHD(){
  phynumber = (gridsetting::dimension>=2 ? 8:7);
  physical_variable_name.push_back("density");
  physical_variable_name.push_back("momentum1");
  physical_variable_name.push_back("momentum2");
  physical_variable_name.push_back("momentum3");
  physical_variable_name.push_back("energy");  
  if(gridsetting::dimension>=2)
    physical_variable_name.push_back("b1");
  physical_variable_name.push_back("b2");
  physical_variable_name.push_back("b3");
  // TODO: transmit needs improvement
  transmit_equ_parameter();
  append_parameter("gamma",5.0/3.0);
  append_parameter("b1constant",0.0);  
}

MHD::~MHD(){
}

void MHD::get_equation_flux(int dim, double* dir, double* var, double* fequ){
  for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
    fequ[tmpi]=0.0;
  if(dim==1)
    get_1d_fx(fequ,var);
  if(dim==2 || dim==3){
    if(std::abs(dir[0])==1.0)
      get_2d_fx(fequ,var);
    if(std::abs(dir[1])==1.0)
      get_2d_fy(fequ,var);
  }
  if(dim==3 && std::abs(dir[2])==1.0)
    get_3d_fz(fequ,var);  
}

double MHD::get_pressure(double* var){
  int dim = gridsetting::dimension;
  double pressure,kinetic_e,magnetic_e;  
  if(dim==1){
    kinetic_e=1.0/2.0*(std::pow(var[1],2)+std::pow(var[2],2)+std::pow(var[3],2))/var[0];
    magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(para["b1constant"],2));
    pressure=(para["gamma"]-1.0)*(var[4]-kinetic_e-magnetic_e);	
  }
  if(dim==2 || dim==3){
    kinetic_e=1.0/2.0*(std::pow(var[1],2)+std::pow(var[2],2)+std::pow(var[3],2))/var[0];
    magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
    pressure=(para["gamma"]-1.0)*(var[4]-kinetic_e-magnetic_e);
  }
  return pressure;
}

void MHD::get_local_vmax(int index){
  int dim = gridsetting::dimension;
  double pressure;
  double var[memory::phynumber];
  for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
    var[tmpi]=memory::variable_tmp[tmpi][index];
  pressure=get_pressure(var);
  if(dim==1)
    memory::local_vmax[index] = 2.0*std::sqrt(para["gamma"]*pressure/var[0]+(pow(para["b1constant"],2)+pow(var[5],2)+pow(var[6],2))/var[0]);
  if(dim==2 || dim==3)
    memory::local_vmax[index] = std::sqrt(para["gamma"]*pressure/var[0]+(pow(var[5],2)+pow(var[6],2)+pow(var[7],2))/var[0]); 
}

void MHD::get_1d_fx(double* fequ, double* var){
  double pressure,magnetic_e, bdotv;
  pressure = get_pressure(var);
  magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(para["b1constant"],2));  
  bdotv=(para["b1constant"]*var[1]+var[5]*var[2]+var[6]*var[3])/var[0];
  fequ[0]=var[1];
  fequ[1]=std::pow(var[1],2)/var[0]+pressure+magnetic_e-std::pow(para["b1constant"],2.0);
  fequ[2]=var[1]*var[2]/var[0]-para["b1constant"]*var[5];
  fequ[3]=var[1]*var[3]/var[0]-para["b1constant"]*var[6];
  fequ[4]=(var[4]+pressure+magnetic_e)*var[1]/var[0]-para["b1constant"]*bdotv;
  fequ[5]=(var[5]*var[1]-para["b1constant"]*var[2])/var[0];
  fequ[6]=(var[6]*var[1]-para["b1constant"]*var[3])/var[0];  
}

void MHD::get_2d_fx(double* fequ, double* var){
  double pressure,magnetic_e, bdotv;
  pressure = get_pressure(var);
  magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
  bdotv=(var[5]*var[1]+var[6]*var[2]+var[7]*var[3])/var[0];
  fequ[0]=var[1];
  fequ[1]=std::pow(var[1],2)/var[0]+pressure+magnetic_e-std::pow(var[5],2.0);
  fequ[2]=var[1]*var[2]/var[0]-var[5]*var[6];
  fequ[3]=var[1]*var[3]/var[0]-var[5]*var[7];    
  fequ[4]=(var[4]+pressure+magnetic_e)*var[1]/var[0]-var[5]*bdotv;
  fequ[5]=0.0;
  fequ[6]=(var[6]*var[1]-var[5]*var[2])/var[0];
  fequ[7]=(var[7]*var[1]-var[5]*var[3])/var[0];
}

void MHD::get_2d_fy(double* fequ, double* var){
  double pressure,magnetic_e, bdotv;
  pressure = get_pressure(var);
  magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
  bdotv=(var[5]*var[1]+var[6]*var[2]+var[7]*var[3])/var[0];    
  fequ[0]=var[2];
  fequ[1]=var[2]*var[1]/var[0]-var[6]*var[5];
  fequ[2]=std::pow(var[2],2)/var[0]+pressure+magnetic_e-std::pow(var[6],2.0);        
  fequ[3]=var[2]*var[3]/var[0]-var[6]*var[7];
  fequ[4]=(var[4]+pressure+magnetic_e)*var[2]/var[0]-var[6]*bdotv;
  fequ[5]=(var[5]*var[2]-var[6]*var[1])/var[0];
  fequ[6]=0.0;    
  fequ[7]=(var[7]*var[2]-var[6]*var[3])/var[0];
}

void MHD::get_3d_fz(double* fequ, double* var){
  double pressure,magnetic_e, bdotv;
  pressure = get_pressure(var);
  magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
  bdotv=(var[5]*var[1]+var[6]*var[2]+var[7]*var[3])/var[0];    
  fequ[0]=var[3];
  fequ[1]=var[3]*var[1]/var[0]-var[7]*var[5];
  fequ[2]=var[3]*var[2]/var[0]-var[7]*var[6];    
  fequ[3]=std::pow(var[3],2)/var[0]+pressure+magnetic_e-std::pow(var[7],2.0);            
  fequ[4]=(var[4]+pressure+magnetic_e)*var[3]/var[0]-var[7]*bdotv;
  fequ[5]=(var[5]*var[3]-var[7]*var[1])/var[0];
  fequ[6]=(var[6]*var[3]-var[7]*var[2])/var[0];
  fequ[7]=0.0;    
}
